Nonlinear dynamics and rheology of active fluids: simulations in two dimensions 
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We report simulations of a continuum model for (apolar, flow aligning) active fluids in two di- 
mensions. Both free and anchored boundary conditions are considered, at parallel confining walls 
that are either static or moving at fixed relative velocity. We focus on extensile materials and find 
that steady shear bands, previously shown to arise ubiquitously in ID for the active nematic phase 
at small (or indeed zero) shear rate, are generally replaced in 2D by more complex flow patterns 
that can be stationary, oscillatory, or apparently chaotic. The consequences of these flow patterns 
for time-averaged steady-state rheology are examined. 

PACS numbers: 87.10.-e, 47.50.-d, 47.63.Gd, 83.60.Fg 
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I. INTRODUCTION 



Active systems in general, and active fluids in particu- 
lar, have recently become a topical research area in soft 
matter physics |lMl9|. An "active particle" is a particle 
which consumes energy from the surrouding environment 
to do work. In many cases this work is used for self- 
propulsion. An active fluid is then a suspension of such 
active particles in an underlying Newtonian fluid. The 
distinguishing feature of activity in this context is the 
ability of the particle to exert forces on the surrounding 
fluid. In the absence of body forces applied externally, 
the simplest perturbation which a particle can impart on 
the fluid is a force dipole. According to the direction of 
the force pair making up the dipole, a single particle can 
be either "extensile" - if the forces are exerted from the 
centre of mass to the fluid - or "contractile" - if it is 
exerted from the fluid to the centre of mass of the parti- 
cle. (For a rodlike particle, extensile forcing ejects fluid 
in both directions along the axis and draws it in around 
the equator; the reverse is true for contractile motion.) 
Also, the dipole can be applied at the centre of drag of 
the object or elsewhere - only in the latter case can this 
active forcing lead to motility (creating a "mover" as op- 
posed to a "shaker" The presence of a force dipole in 
active fluids defines a director at the particle level even 
in cases where the organism is not rodlike. Thus one 
can introduce a natural order parameter which charac- 
terises the magnitude of local orientational order in the 
fluid. As a result, in concentrated active fluids one gener- 
ically expects an isotropic-nematic transition, and this is 
why such systems are commonly described by means of 



equations of motion that closely resemble those of liquid 
crystalline materials. 

A suspension of bacterial swimmers, such as B. sub- 
tilis, is an example of an extensile active fluid, whereas 
a suspension of algae like Chlamydomonas is contractile. 
At a quite different length scale, the same equations can 
be used to describe an actomyosin solution containing fil- 
aments and motors. In this case the active motion of the 
motors, which resemble moving cross-links on the net- 
work of filaments, creates a contractile effect. Such acto- 
myosin networks represent a simplified picture of the cy- 
toskeletal structures underlying the motility of eukaroytic 
cells [sl, 3 . In this work we prefer the term "active fluids" 
to "active gels" for all the systems under study, although 
the latter term is widely used. (In fact not all these ma- 
terials have a markedly viscoelastic rheological response, 
as the term "gel" suggests; we shall discuss this later 
on.) The very definition of active particles implies that 
they are far from equilibrium; unsurprisingly therefore, 
nonequilibrium phenomena are a dominant theme in the 
physics of active fluids. 



For instance, it has been predicted theoretically 
that active fluids should undergo a nonequilibrium phase 
transition between a "passive" quiescent phase, where 
the motion of each of the particles are basically uncorre- 
lated and the coarse grained mean velocity field is uni- 
form and zero, as in conventional passive unforced fluids, 
and an "active" phase, in which long-range correlations 
lead to a non-zero "spontaneous flow" in steady state. 
This transition has also been seen experimentally: there 
the spontaneous flow patterns have been mapped out by 
following the dynamics of individual bacteria in a thin B. 
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subtilis film, where the local concentration of bacteria was 
very high [6|. What triggers the transition between the 
quiescent and the "active" phase in these bacterial fluids 
is simply the increase in concentration - which controls, 
among other things, the magnitude of an "activity" pa- 
rameter which appears in the equations of motion that we 
will consider later on. The bacterial flow patterns invari- 
ably were seen to involve vortices and swirls of bacteria, 
similar to those which are observed in aerobic bacteria 
which move around to get into contact with oxygen 0. 
These large-scale motions are sometimes referred to as 
"bacterial turbulence" because the flow field resembles 
that of a Newtonian fluid at high Reynolds number, in 
the turbulent regime. However, it should be kept in mind 
that, unlike the turbulent flow of Newtonian fluids, that 
of active fluids occurs at essentially zero Reynolds num- 
ber, so that inertia plays no role. 

Although less characterised experimentally as yet, the 
rheology of active fluids is also expected to display a very 



intriguing phenomenology [SMlOL Il5| . Extensile and con- 
tractile fluids lead to very different rheological responses. 
It was predicted and later on confirmed by simulations 
in ID that contractile active gels should show an almost 
solid-like behaviour when sheared, and that if left free 
to reorient at the boundary they should possess a non- 
zero yield stress Extensile fluids on the other 
hand should have in ID a "negative yield stress" causing 
a window of apparently superfluid rheology (in which a 
nonzero macroscopic shear rate can be accommodated at 
zero stress) close to the transition between the isotropic 
and the nematic phase. 

An outstanding problem in the continuum theory of 
active fluids is that most of the calculations done so far 
assume a ID (or rather a quasi- ID) geometry, in which 
the variations of the orientational order parameter and 
the flow field are limited to just one dimension. (The ori- 
entational order itself lives in a 2D or 3D space.) While 
this is enough to describe the spontaneous flow transition 
and to map out rheology curves, it is legitimate to ask 
how much the ID predictions are confirmed by fully 2D 
or 3D simulations. (The same would hold in most fluid 
dynamics problems involving the occurrence of hydrody- 
namic instabilities; spontaneous flow is just one of these.) 
A small number of simulations in 2D for extensile spon- 
taneously flowing fluids, within an active but isotropic 
phase, have already shown that flow patterns can indeed 
be significantly different in 2D than in ID. For instance, 
while in ID a spontaneous net flow is the first state found 



on entering the active phase, in 2D rolls and turbulent 
flow were instead observed, which are more closely re- 
lated to the experiments in Ref. [gI, Q . 

Although, as stated above, there are few works in 2D 
simulating active fluids using the hydrodynamic contin- 
uum equations of [sj as we pursue here, several other 
simulation avenues have been pursued in 2D. These in- 
clude a two-phase model, developed by Wolgemuth [21| in 
a context specific to bacteria; this focuses on the coupling 
between activity and number density, which we neglect in 
this paper. Other contributions have involved simulation 
of discrete swimmers in two or three dimensions |17l. l22- 
25]; however, such approaches become increasingly dif- 
ficult as one approaches the case of dense swarms for 
which the hydrodynamic description (without coupling 
to number density) was primarily developed. Several of 
these methods have shown onset of bacterial turbulence 
or development of coherent structures of various kinds, 
but it is not yet clear, for instance, which of these require 
the coupling of activity to number density. 

Our aim in this paper is therefore to study in more de- 
tail the simplest continuum models for 2D flowing states 
of active fluids, focusing on the extensile case, close to or 
within the nematic phase, both in the absence of any forc- 
ing and when subjected to a small shear rate. Besides, 
we evaluate the role of boundary conditions by comparing 
simulations in which the nematic order parameter field 
is fixed on the surface (to ensure planar alignment of the 
associated director field) to others in which the director 
is free to rotate at the boundary planes. We confirm that 
2D simulations of spontaneous flow patterns in the active 
phase markedly differ from ID patterns. We also show 
that the rheology curves may differ when the ID approx- 
imation is lifted. In essence, this paper extends to 2D the 

work presented for ID systems in our earlier paper with 
n 

Orlandini and Yeomans [9]. Naturally one can ask what 
would happen if one generalizes further to the fully 3D 
case. A comprehensive study of this case would be com- 
putationally exhausting, and we defer it to future work. 
However, some initial progress has recently been made 
using stability analyses to guide a selective exploration 
of parameter space |2Ql] . 

The rest of this paper is structured as follows. In Sec.lTTl 
we review the equations of motion of apolar active flu- 
ids, to which we restrict ourselves here. In Sec. IIII Al 
we specify for our simulations the geometry and bound- 
ary conditions; the parameter values and units used; and 
(briefly) the algorithmic methods involved. In Sec. IIVI 
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we summarise and extend slightly our earlier results 
for the ID hydrodynamics and rheology of these active 
suspensions. The results of our new 2D simulation study 
are then presented in Sec. |V]- note that we focus on ac- 
tive extensile materials which lead to spontaneous flow in 
a quasi- ID geometry, for rod-like active particles 
Sec. |Vl] contains conclusions and perspectives for future 
study. 

II. MODEL EQUATIONS 



Following standard procedures [IQ^ we first employ a 
Landau - de Gennes free energy J- = J fdV to describe 
an equilibrium isotropic to nematic transition in a ma- 
terial without activity. The free energy density can be 
written as a sum of two terms, / = /i + /2- The first is 
a bulk contribution 



(QL)'(i) 



while the second is a distortion term, which we take in a 
(standard) one-constant approximation as ^[ 



/2 = f (d^C 



'a(3 



(2) 



In the equations above, Aq is a constant measuring the 
relative contributions of the bulk and the distortion term, 
7 controls the magnitude of order (it may be viewed as an 
effective temperature or concentration for thermotropic 
and lyotropic liquid crystals respectively), while K is an 
elastic constant. The resulting free energy density f i s 
standard to describe passive nematic liquid crystals [26| . 
Here and in what follows Greek indices denote cartesian 
components and summation over repeated indices is im- 
plied. 



The equation of motion for Q is taken to be \27V\2 



{dt^u- V)Q - S(Vu, Q) = FH + AQ 



(3) 



where F is a collective rotational diffusion constant. We 
have added a term proportional to A, which is one of two 
"activity parameters" that become nonzero when activity 
is switched on. However, this term is easily absorbed into 
a redefinition of the coefficient 7 in Eq. ([Tj) and from 
now on we suppress it: A = 0. The form of Eq. was 
suggested on the basis of symmetry in Refs. |3| and 
derived starting from an underlying microscopic model 
in Ref. [ij]. The first term on the left-hand side of 
Eq. (|3]) is the material derivative describing the usual 
time dependence of a quantity advected by a fluid with 



velocity u. This is generalized for rod-like molecules by 
a second term 



S(Vu,Q) = (eD + cj)(Q + I/3) 
+ (Q + I/3)(eD-cj) 
- 2^(Q + I/3)Tr(QVu) 



(4) 



where Tr denotes the tensorial trace, while D = (Vu + 
Vu^)/2 and uj = (Vu — Vu^)/2 are the symmetric part 
and the anti-symmetric part respectively of the veloc- 
ity gradient tensor Vu^js = d^u^- The constant de- 
pends on the molecular details of a given particle and 
controls whether the passive material is flow-aligning or 
flow-tumbling. To isolate the nolinear effects of activity 
we choose flow-aligning materials in this paper (bearing 
in mind that flow-tumbling ones often show unsteady flow 
behaviour even without activity). The first term on the 
right-hand side of Eq. ^ describes the relaxation of the 
order parameter towards the minimum of the free energy. 
The molecular field H which provides the force for this 
motion is given by 



(5) 



The fluid velocity, u, obeys the continuity equation for 
an effectively incompressible fluid, V.u = 0, and also the 
corresponding Navier-Stokes equation, 

p{dt + Uf3df3)Ua = d^{Iia(3) + r]9/3{daU/3 + 9/3^/^) (6) 

Here p is the fluid density, 7^ is a viscosity and 11^(3 = 



-rrpassive 



H^^^^®. Note that this viscosity rj in Eql6] 



is isotropic, and that of the solvent in which our ac- 
tive particles are suspended. This does not preclude the 
emergence of an anisotropic macroscopic viscosity as a re- 
sult of the coupling to those particles. The stress tensor 
H^^^^^^^ necessary to describe ordinary LC hydrodynam- 
ics without activity is: 



T-i-passive 

^^a(3 — 



Po^a(3 + 2^((5q,/3 + -Sa(3)Q-feH^e (7) 



5F 



The nematic free energy gives an effectively constant con- 
tribution to the fluid pressure Pq (whose final spatial vari- 
ation is determined by the fluid incompressibility condi- 
tion). The active stress, which unlike the passive one 
cannot be derived from any underlying free energy func- 
tional, is given to leading order by 



n active /-/o 
a/3 = -CVa/S 



(8) 
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where is a second activity constant [l|, [isj . Note that 
with the sign convention chosen here C > corresponds 
to extensile rods and ^ < to contractile ones [l|; in 
either case this term does not merely renormalize the 
equations for a passive liquid crystal but fundamentally 
alters their form. Accordingly it is a key control param- 
eter in the continuum description of active nematics. As 
in Eq. [3l the explicit form of the active contribution to 
the stress tensor entering Eq. [6] was proposed on the 
basis of a symmetry analysis of a fluid of contractile or 
extensile dipolar objects in It was also derived by 
coarse graining a more microscopic model for a solution 
of actin fibers and myosins in Ref. [14 1. 

The equations of motion chosen above address the case 
of active nematics, by which we mean particles whose 
local ordering is apolar. This means that locally one 
has a preferred orientation of the force dipole but not of 
any vector field such as the mean swimming direction of 
motile particles. The equations of motion in the latter 
case are yet more complicated (see e.g. [lol) but because 
the active stress takes the same form as above, for rhe- 
ological purposes are expected to yield broadly similar 
results. In this work, for simplicity, we study solely the 
apolar case. 



III. SIMULATION DETAILS 

In this Section we describe the geometry and boundary 
conditions used for our numerical work and also discuss 
the units and parameter values chosen, and the issue of 
numerical convergence. Note that (primarily for histori- 
cal reasons) we used finite difference (FD) for free bound- 
aries and lattice Boltzmann (LB) for fixed ones; this is 
considered further in Appendix A. 



A. Geometry and boundary conditions 

We consider a slab of active fiuid sandwiched between 
fiat parallel plates located at y = O^Ly. The fiuid ve- 
locity and order parameter tensor can vary in both x 
and while in the third dimension z we assume transla- 
tional invar iance. Each plate has length Lx and periodic 
boundary conditions are imposed in the x direction. In 
general the upper plate is taken to move in the positive 
X direction at a constant speed jLy^ although many of 
the results presented below are for the case with no ex- 
ternally applied shear, 7 = 0. At the plates we impose 



boundary conditions of no slip and no permeation for the 
fiuid velocity. For the molecular order parameter, we will 
study two different boundary conditions: 

• So-called "free" boundary conditions, in which the 
order parameter tensor at the wall can take any 
value but must satisfy a zero-gradient condition in 
the direction normal to the wall: 



dyQa(3 = 0, at y = {0,Ly} y a, (3 



(9) 



Free boundaries are believed to give the fastest conver- 
gence to bulk behavior by minimizing the effects of the 
confining walls on the order parameter dynamics. 

• So-called "fixed" boundary conditions, in which the 
director field is anchored along the x direction par- 
allel to the wall, and the order parameter tensor 
Qai3 has the equilibrium form for a uniaxial state 
with this director: 



Qa/3 q riaUis — 
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where 



(10) 



(11) 



is the magnitude of order in the surface, which we 
take to be equal to the one in the bulk. 

The choice of fixed boundary conditions is motivated 
by the behavior of non-active liquid crystals which fre- 
quently have specific anchoring interactions that lock in 
the director field at the boundary. Their appropriateness 
for active nematics is less well established (particularly 
for bacterial suspensions, although some sort of anchor- 
ing remains plausible for actomysin gels), which is why 
a comparison of the two types of boundary condition is 
appropriate here. In fact, as we will see below, our main 
conclusions concerning the fiow patterns and rheology are 
relatively robust to the choice of boundary conditions (al- 
though these certainly infiuence some of the details). 



B. Units and parameter values 

Throughout our study we use units of length in which 
the gap between the plates Ly = 1; units of time in 
which the model's underlying microscopic timescale r = 
l/{AoT) = 1; and units of mass such that the free energy 



parameter Aq, which has the diniensions of a stress, hke- 
wise obeys Aq = 1. Ah runs use a value of the I/N control 
parameter 7 = 3.0 which lies within the nematic phase 
for a system without activity; as stated previously this 
value effectively absorbs the activity parameter A, which 
we have set to zero. Throughout we set ^ = 0.7, which 
corresponds to a typical value for flow- aligning molecular 
nematics [26]; such nematics generally comprise rod-like 
molecules of modest aspect ratio. In any case we ex- 
pect that, so long as one does not leave the flow- aligning 
regime, our results should be relatively insensitive to this 
choice. In the system of units just defined, we choose 
for the Newtonian viscosity a value 77 = 0.567, which is 
much larger than the viscoelastic one without activity, 
hence leading to a conventional Newtonian rheology in 
the passive phase. 

With the above choices, and denoting by 7 the mean 
shear rate imposed by the relative motion of the two 
plates, the parameter values that remain to be varied 
between runs are the activity level and a microscopic 
lengthscale on which elastic distortions compete with 
bulk free energies. This obeys / = ^/K/Ao = (where 
the latter equality holds in our system of units only). Ac- 
cordingly, we will quote values for the parameters (C, I) 
in each figure caption of our results sections below. A 
third parameter, the Reynolds number (Re) which con- 
trols the relative strength of inertial to viscous terms in 
the Navier-Stokes equation, is always small or zero (see 
below) . 

C. Numerical convergence 

In a class of systems whose flow frequently is (or ap- 
pears to be) chaotic, one cannot expect genuine conver- 
gence of numerically calculated velocity patterns with 
respect to the time-step At and mesh scale (Ay^Ax). 
We strive to ensure our results are converged, in the 
sense that further refinement gives no significant change 
to the type of behavior observed (except very close 
to parameter values at which there is transition from 
one regime to another), nor to macroscopic quanti- 
ties like the time-averaged stress. In the units cho- 
sen above, we found this to typically require values of 
{Ay, Ax, At) = (1/100, 1/100, 0.34) in our LB simula- 
tions and (A^, Ax, At) = (1/128, 4/512, 0.05) - or some- 
times {Ay, Ax, At) = (1/256, 4/1024, 0.05) - in our FD 
simulations. For sufficiently small /, one expects the 
boundary conditions to be unimportant so that (mod- 




FIG. 1: Homogeneous constitutive curves. 

ulo the minor technical differences discussed above) the 
results from FD and LB should approach one another in 
this limit. 

IV. RESULTS IN ZERO AND ONE DIMENSION 

In this section we briefly recap some results of Ref. [sl 
for the hydrodynamics and rheology of an extensile fluid 
in less than two dimensions. We discuss first the case of 
a homogeneous (OD) imposed shear flow, and then recall 
the results of calculations that allow spatial variations 
in one spatial dimension (ID), y, assuming translational 
invar iance in both z and x. 

The homogeneous constitutive curves are shown in 
Fig. [1] for three different values of positive (extensile) ac- 
tivity, C > 0- As discussed in Ref. 0, the vertical drop 
at the origin arises from the alignment by weak flow of 
the nematic director at the Leslie angle (this occurs, for 
a passive system, throughout the nematic phase). This 
alignment produces an active stress tensor whose xy com- 
ponent for extensile particles is of opposite sign with re- 
spect to shear rate. This negative stress contribution is 
proportional to ( and depends on the sign of the imposed 
flow 7 but not its magnitude (when this is small). This 
discontinuous variation creates a downward step in the 
flow curve, which is expected to allow the spontaneous 
formation of shear bands, even when no external shear is 
applied at the system's boundaries, 7 = |9|. 

Accordingly, we perform a series of unsheared ID runs 
in which spatial variations are allowed in y (only). Each 
is initialised with the tensor Qa(3{y) having a uniform {y- 
independent) degree of uniaxial order, and its long axis in 
the xy plane at an angle = 22° tanh((0.5 — y)/l) to the 
X direction. This favours the formation of just two shear 
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FIG. 2: Top: shear banded (or quiescent) profiles for ( = 
0.01. Values of / =0.0005, 0.00071, 0.001, 0.00113, 0.0014, 
0.002, 0.0028, 0.004, 0.00476, 0.0057, 0.00673, 0.0073, 0.008, 
0.00872, 0.00951, 0.016 increasing for decreasing j{y = 0). 
Inset: shear rate = 0) at the left edge of the cell extracted 
from such profiles, for various values of shown in the 
master scaling representation discussed in the main text. 
Bottom: Phase diagram for ID runs showing the region of 
spontaneous shear banded flow (closed circles) and the region 
in which an initially heterogeneous state decays back to a 
quiescent state of zero flow (open circles). The dashed line is 
the power law /* = aC^^^ using the value of the intercept a 
extracted from the master scaling plot in the inset to the top 
figure. 

bands. (For random initial conditions, multiple bands 
can in principle form in this planar flow geometry.) The 
code was then evolved to steady state. 

In the limit / ^ for any value of ( we find the steady 
state to comprise two coexisting bands of equal and op- 
posite shear rates ±7*(("), these being the values at which 
the two positively sloping branches of the homogeneous 



constitutive curve intercept the 7 axis (see Fig. [2j). The 
bands are separated from each other by a slightly diffuse 
interface of thickness I. Increasing / in a series of runs 
at any fixed value of ( eventually eliminates this state 
of spontaneous flow, with quiescence being restored at a 
critical value /* = a^^/^ where a ~ 0.885. This is seen in 
the master scaling curve of j{y = 0,t ^ oo)/7* versus 
IjC}^^ in the inset of Fig. [2l top, and by the dashed line 
separating quiescence (open circles) from spontaneous 
banded states (filled circles) in Fig. [21 bottom. 

The results just presented were found with free bound- 
ary conditions. Simulations of the fixed boundary case 
give a more complex behavior, particularly at the larger 
values of /, where the anchoring condition at the wall can 
compete with the elastic director distortions required to 
maintain the shear-banded state. For a fuller discussion 
of such effects in ID see 

V. 2D RESULTS: EXTENSILE SYSTEMS 

The remainder of the paper concerns 2D simulations 
that allow spatial variations parallel to the plates (along 
x) as well as perpendicular to them (along y). Further- 
more we focus on extensile fluids, which lead to spon- 
taneous flow in ID (contractile systems will be treated 
elsewhere). 

A. Unsheared systems 

We discuss first systems that are not subject to any 
externally applied shear flow, treating the free and fixed 
boundary condition cases in turn. (Systems with applied 
shear will be addressed in Sec. IV Bl below.) 

1. Free boundary conditions 

Using the free boundary condition, Eqn. [9] above, we 
performed a series of 2D runs at the values of ^ shown 
by squares in Fig. [3l In each run we input as an initial 
condition the final state of the corresponding ID run of 
Sec. IIVI above. (All these ID runs had themselves used 
free boundary conditions.) The dashed line in Fig. [3] is 
copied directly from Fig. [21 bottom. Accordingly, all runs 
above this line had a quiescent initial condition. All those 
below it started with coexisting shear bands of equal and 
opposite shear rate. To this ID intial condition, we also 
added a 2D random component (the noise distribution 
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FIG. 3: Phase diagram for 2D runs with free boundary condi- 
tions, each denoted by a square. Empty squares: quiescence. 
Elsewhere initial shear banded profile gives way to a state 
of 2D domains that is: steady (crosses); oscillatory (shaded 
square); aperiodic (filled squares). 
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FIG. 4: Throughput versus time for / = 0.002 and C = 0.004 
(left), C = 0.04 (right). 

was uniform, and could be either positive or negative for 
each component of the order parameter tensor) of tiny 
ampHtude. This initiaHsation procedure allowed a study 
of the linear regime kinetics for the initial growth (or 
decay) of 2D perturbations at early times, as well as the 
ultimate attractor that is attained at long times. 

For all runs located above the dashed line in Fig. [3l the 
initial 2D perturbation decayed to zero as a function of 
time, showing the homogeneous quiescent state to be lin- 
early stable. In all runs below the line, for which the ID 
base state was shear banded, the initial 2D perturbation 
grew in time. 

In each case, this destabilisation of the ID initial state 
resulted at long times in a more complicated state with 
2D structure. A macroscopic signature of this evolution 
is the decay of the gap- averaged throughput, which was 
non-zero in the initial V-shaped ID velocity profile, to 



zero (on average) in the final 2D state; see Fig. HI (Note 
that by convention we choose the V-shaped velocity pro- 
file in the initial state to correspond to a positive through- 
put.) 

At long times, the signal of integrated amplitude 
versus time (for all qx modes in any given run) set- 
tled either to steady state, or to an oscillatory attrac- 
tor, or to an aperiodic, apparently chaotic, attractor. 
The word "apparently" is used here because we do not 
measure Lyaponov exponents and therefore do not dis- 
tinguish true chaos from quasiperiodic motion. (From 
now on, though, we neglect this distinction and use the 
terms 'chaotic' and 'aperiodic' interchangeably.) We dis- 
tinguish these three different dynamical regimes by the 
type of filling of the squares in Fig. [3] (crosses=steady; 
shaded=oscillatory; solid- filling=chaotic.) Representa- 
tive snapshots of the full 2D state at a long time after the 
start of each run, once the system has attained this ul- 
timate attractor, are shown in Figs. [5] and [6] respectively 
for runs performed along the horizontal and vertical thin 
dotted hues in Fig. [H As can be seen, steady roll-like 
states tend to dominate towards the upper left of the 
unstable regime in the C, / plane, giving way to chaotic, 
turbulent-like states at the lower right. 

It is remarkable that the ID shear bands are imme- 
diately unstable towards a static "roll-like" flow pattern 
with variation in both x and y. Given this, however, the 
general progression from steady rolls via oscillations to 
chaotic flow on increasing activity at fixed / (that is, fixed 
ratio of sample dimension to microscopic length, L/l) is 
fairly natural - similar behavior was reported anecdo- 
tally for one parameter set in |31|. Notably, the same 
progression is seen at fixed by instead decreasing / - 
equivalent to increasing L//, the ratio of the sample size 
to the molecular length scale. Thus not only is the in- 
stability of the quiescent state delayed in small systems 
(as it is in ID, [9|), but also the subsequent transitions to 
oscillation and chaos are likewise delayed. This reflects 
the high energy cost of creating inhomogeneous director 
fields in systems that are not many times larger than the 
elasticity length /; this cost stabilizes the quiescent state 
but can be overcome, for all the / values studied here, by 
increasing the activity sufficiently. 

2. Fixed boundary conditions 

We now turn to the case of fixed boundary conditions. 
FiglJl shows the phase diagram in this case, plotted in the 



m 



\W4 



m 




m 



eiliiii! M 

4, 




> wm 



FIG. 5: Snapshot at long time of 2D runs for C = 0.002,0.004,0.01,0.02,0.04,0.1 and / = 0.002 (free boundary conditions, 
no external shear throughout). Greyscale on left shows Qxx] while arrows on right show the fluid velocity. The x direction is 
horizontal, y vertical. 



same way as Figl3]for free BCs. Just as we found there, 
the simple ID banded state (creating a V-shaped velocity 
profile) is never seen. As expected, in the regime of small 
/ (large L/l) the BCs become relatively unimportant and 
the behavior as a function of ^ very similar to that with 
free BCs. However, at larger / (i.e., smaller system sizes 
in units of the microscopic length), the choice of BCs 
plays a larger role. 

In ID the fixed BCs, by constraining the director to 
lie in the flow plane, somewhat suppress the instability 
towards the banded state relative to free BCs [9|. How- 
ever, the ID banded state is again unstable; moreover, 
the presence of additional instability modes pushes the 
boundary of the passive phase back towards lower activ- 
ity. 

Specifically, to numerical accuracy we discern the fol- 
lowing progression at fixed / = 0.002. If the activ- 



ity is smaller than, but close to, the critical threshold 
in ID, the system forms steady convection rolls (seen 
previously in [3l|). Fig. [8^ gives a steady-state snap- 
shot of Qxx and of the fluid flow in this regime. This 
steady state is achieved in about 200,000 LB timesteps 
which corresponds to time 6.755 x 10^ in our units (i.e., 
t = 6.755 X lO^r). On the other hand, if the activity 
is raised to just beyond the level needed to create spon- 
taneous flow in ID, then the rolls break and form what 
looks like one or more tilted bands (Fig. [SJd). These 
tilted bands are again stationary. Additional simulations 
in a box with a 4:1 aspect ratio (not shown) suggest that 
these may be stabilised by the reduced geometry and box 
size used here - consistently with the fact that we do not 
observe these states with free boundary conditions in a 
4:1 box (Section VAl). 

If the activity increases further (Fig. [8t, C = 0.008) a 






FIG. 6: Snapshot at long time of 2D runs for C = 0.1 and, from top to bottom, / = 0.016,0.008,0.004,0.002 (free boundary 
conditions, no external shear throughout). Greyscale on left shows Qxx] while arrows on right show the fluid velocity. The x 
direction is horizontal, y vertical. 



roll-like state that forms initially breaks into an appar- 
ently undulating band which then oscillates. This looks 
quite similar to Fig. [6] (second panel from top) which is 
however already chaotic (perhaps because of the larger 
aspect ratio). A further increase to = 0.01 appears to 
stabilise multiple bands with small vortices in between 
(Fig. [8]i). These states are quasi- ID in nature, in the 
sense that the main variation occurs along the velocity 
gradient direction, yet were not seen in our strictly ID 
simulations. For yet larger activity, we have an appar- 
ently chaotic flowing state, similar to that reported above 
for free BCs (FiglH^), also seen previously in Ref. [3l|. 



It should be noted that the threshold for spontaneous 
flow in ID and 2D are different with fixed boundary 
conditions. This is a difference with respect to the free 
boundary condition results discussed before, which is due 
to boundary stabilisation of the passive phase. For rela- 
tively large values of //L, the boundaries provide layers 
within which the ordering is fixed by the anchoring, and 
this provides the observed stabilisation. As expected, we 
find that as the value oil/ L decreases, so does the differ- 
ence between the thresholds in ID and 2D (see Fig. [7]). 



B. Sheared systems 

Here we address systems that are subject to a shear 
flow applied externally at the cell boundaries. As before 
we consider in turn the cases of free and fixed boundary 
conditions. 



1. Free boundary conditions 

For very small 7, one expects the effect of flow to be 
perturbative on the underlying zero-shear states found 
above. In a strictly ID flow showing shear bands, the 
effect is merely to shift the band interface, creating a 
nonzero mean shear rate without altering any structure 
within the bands. This degeneracy allows the system 
to accommodate a small macroscopic shear flow without 
developing a stress (called "superfluidity" in [9]). In ID 
we also found that increasing the shear rate beyond the 
superfluidity window caused elimination of one of the two 
shear bands, thereby restoring a homogeneous shear flow. 

In 2D, with a more complicated inhomogeneous flow 
in the absence of bulk shear, it is far from clear whether 
a similar degeneracy ought to exist. If it does not, one 
expects a finite shear rate to be accompanied by a finite 
stress, so that the superfluid region is replaced by one 
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FIG. 7: Top: Phase diagram for ID runs with fixed boundary 
conditions, each denoted by a circle. Empty circles denote 
quiescence, filled ones spontaneously flowing states. Bottom: 
Phase diagram for 2D runs with fixed boundary conditions, 
each denoted by a square. Empty squares: quiescence. Else- 
where the stable ID profile with a small inhomogeneity gives 
way to a state of 2D domains that is: steady (crosses); oscil- 
latory (shaded square); chaotic (filled squares). 



of finite (zero-shear) viscosity. Nonetheless one might 
expect restoration of homogeneous flow at high enough 7, 
either by a sharp transition or only gradually as 7 ^ 00. 

To address the behavior in 2D we use for reference the 
unsheared phase diagram in the ((", I) plane of Fig. [3l 
focusing on three locations in it: / = 0.002 and ( = 
0.001,0.005,0.04, for which the unsheared state com- 
prised simple rolls, wavy rolls and turbulence respec- 
tively. (Recall Fig. [5l) At each of these locations we 
perform a series of runs for increasing values of the ap- 
plied shear rate 7. At each shear rate separately we 
perform shear startup from a state of no flow, with 
an initial condition in which the order parameter ten- 
sor is taken to be everywhere uniaxial for simplicity, 
Qa/3 = \Si{nocnf3 — (5^/3/3), and with its principal axis 
assumed to reside in the xy plane. To avoid biasing the 
system into any particular final state, at each individual 



grid point separately we assign a random value for the 
order parameter (drawn from a flat top-hat distribu- 
tion); and for the angle of the principal axis (from a fiat 
distribution between and 27r). 

Our results are shown in Figs. [9] to [Til For low ac- 
tivity C = 0.001, the static roll pattern observed with- 
out shear is destroyed even at the smallest shear rates 
studied, whereas at high enough shear rates a homoge- 
neous flow is recovered and the shear stress reverts to 
the one calculated in ID. The flow pattern at very low 
shear rates is unsteady and is spatially nonuniform with 
no evident periodicity in either space or time. The shear 
stress takes both positive and negative values. The un- 
steadiness makes it very difficult to determine a reliable 
value for the average shear stress at small 7 (error bars 
in Fig. [9] are standard deviations of the mean from the 
time series). 

At C = 0.005, the 'wavy roll' state observed in the 
case of zero bulk shear (closely resembling that of Fig. [5]) 
is again broken up at the smallest shear rates studied. 
This is not surprising since this state has some form of 
(defect-ridden) layerwise order along the flow direction 
which is inconsistent with the imposed shear geometry. 
In contrast to the preceding case, the system now shows 
clear evidence of a finite viscosity at low shear rates (with 
much smaller stress fluctuations at low 7). Again, the 
system approaches a homogeneous laminar flow at high 
enough flow rates where the stress similarly attains the 
value predicted from the ID calculation in this regime. 

Finally, for = 0.04, the chaotic state present initially 
is perturbed only marginally at small strain rates. The 
data is consistent with a finite effective viscosity at low 
strain rates, but unlike the previous two cases the flow 
curve has noticable upward curvature here (so a ^ 
with p > 1 cannot be ruled out). At high strain rates the 
flow becomes progressively more homogeneous and the 
shear stress intersects the ID curve, signifying a return 
to a laminar ID flow. 



2. Fixed boundary conditions 

We now turn to the case of fixed boundary conditions. 
In ID the results are very similar to the free bound- 
ary case [sl. The main difference between fixed and free 
boundary conditions in ID is that the superfiuid window 
is somewhat narrowed by the constraint on the order pa- 
rameter (which inhibits the formation of a band interface 
close to a wall) and can disappear altogether for large 
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FIG. 8: Snapshots corresponding to simulations with extensile fluids and fixed boundary conditions, for a system with Lx — 
Ly — 1 and periodic boundary conditions along the flow direction. The snapshots show grayscale plots of Qxx (left column), 
and the velocity profile (right column). The value of the activity is (from top to bottom): 0.002 (a), 0.004 (b), 0.008 (c), 0.01 
(d), 0.08 (e); throughout, / = 0.002. 
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FIG. 9: Top: Flow curve for C = 0.001,/ = 0.002. Lower: 
representative state snapshot at a late time (on the final at- 
tractor) for each of the lowest seven shear rates. (Shear rate 
increasing in snapshots downwards. The x direction is hori- 
zontal, y vertical.) The solid line shows the OD homogeneous 
constitutive curve. 

enough value of / 0. 

As done in the preceding Section focussing on free 
boundary conditions, we here report flow curves and typ- 
ical stress snapshots found for different values of activity, 
which we chose as representative of the steady rolls and 
the turbulent regimes in the unsheared phase diagram. 




FIG. 10: Top: Flow curve for C = 0.005, / = 0.002. Lower: 
representative state snapshot at a late time (on the final at- 
tractor) for each of the lowest seven shear rates. (Shear rate 
increasing in snapshots downwards. The x direction is hori- 
zontal, y vertical.) The solid line shows the OD homogeneous 
constitutive curve. 



The flow curve observed when shearing the static roll 
pattern is shown in Fig. [121 One important qualitative 
difference with respect to the free boundary case is that 
there appears to be a linear regime in this case, in which 
the viscosity of the rolls is well defined, and larger than 
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FIG. 11: Top: Flow curve for C = 0.04,/ = 0.002. Lower: 
representative state snapshot at a late time (on the final at- 
tract or) for each of the lowest seven shear rates. (Shear rate 
increasing snapshots downwards. The x direction is horizon- 
tal, y vertical.) The solid line shows the OD homogeneous 
constitutive curve. 



77 - see Fig. [131 The extent of this linear regime is very 
small, and quite likely it shrinks with increasing Lx - 
which may explain this qualitative difference between the 
free and fixed boundary observations. On the other hand, 
in good agreement with the free boundary case, at high 
enough shear rates a ID flow curve is recovered. At low 



7, again as with free boundaries, the shear stress takes 
both positive and negative values. 

The behaviour shown in Fig.[ll]for ( = 0.04, is selected 
as typical for flow curves starting from the unsheared 
turbulent regime. The flow curve is quantitatively in 
good agreement with its counterpart for free boundary 
conditions shown before. The data at small shear rate 
are consistent with a small yield stress in this case - 
this is reasonable as turbulent flows may be expected to 
dissipate more than rolls, hence have a larger viscosity. 
The linear regime has also disappeared in the turbulent 
regime. 



VI. CONCLUSION 

To conclude, we have presented numerical calculations, 
by means of both finite difference and Lattice Boltzmann 
simulations, of the dynamics and flow behaviour of active 
extensile fluids, such as, for instance, concentrated bacte- 
rial suspensions. After checking consistency between the 
two methods we have used them to explore the differ- 
ent active flows associated with free and fixed boundary 
conditions. 

In the unsheared case, we found that the transition to 
spontaneous flow occurs significantly differently than in 
ID, consistently with previously reported 2D simulations 
performed with selected parameters [l9|. Instead of ac- 
tive bands, we find that the spontaneously flowing phase 
typically consists of rolls or turbulent flow. In some runs, 
we were nevertheless able to stabilise tilted active bands, 
which are however different in shape from the purely ID 
ones. 

The 2D flow curves are also rather different from the 
ID ones. Our finite difference simulations with free 
boundary conditions were performed with a 4:1 aspect 
ratio between the flow and the flow gradient directions, 
and this leads to a vanishing or very small linear regime. 
With fixed boundary conditions, on the other hand, we 
find that rolls have a well-defined viscosity in the linear 
regime, and this is larger than the nominal Newtonian 
viscosity of the fluid. This is in agreement with the ex- 
pectation that 2D vortices lead to enhanced dissipation 
with respect to ID bands. 

Fully 3D simulations of this problem would undoubt- 
edly be worthwhile, but present serious computational 
challenges \19\ . An intermediate step to a full 3D simula- 
tion is to maintain a two dimensional simulation domain 
but allow the simulated velocities and director fields to 
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FIG. 12: Top: Flow curve for C = 0.001,/ = 0.002. The 
dashed curve is part of the homogeneous constitutive curve. 
Lower two panels: representative state snapshot at a late time 
for 7 = 10~^ (second row) and 7 = 10"^ (third row . 

be fully 3D. This is actually what we do in our LB stud- 
ies (see Appendix) although in all the states discussed 
so far the flow and director field remains two dimen- 
sional, provided that the initial configuration is purely 
two-dimensional. If this is not the case, a 3D flow may 
result, where both the director and the velocity field ac- 
quire an out-of-plane component (an example is shown 
in Fig. [T5]). Intriguingly this has a linear viscosity which 
is smaller than the Newtonian one, in line with the ID 
results and dissimilar to the main trends observed in 2D. 
The details of this out-of-plane flow mode will be pursued 
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FIG. 13: Effective viscosity (normalised by 77) versus time for 
an active gel with / = 0.002, ( = 0.001 with fixed boundary 
conditions. The solid line refers to 7 = 0.0002, the dashed one 
to 7 = 0.001. The steady state effective viscosities are very 
similar in both cases, proving that rolls have a linear rheology 
regime. 

elsewhere. 

All our findings can be summed up to say that the dy- 
namics and rheology of active extensile fluids is extremely 
nonlinear, and that one should take any prediction com- 
ing from purely ID flow, which are typically considered 
in analytical calculations, with care, as 2D results are 
often different and more complicated. Nevertheless, our 
results offer firm predictions for macrorheological exper- 
iments on thin 2D samples, which could be performed 
by e.g. shearing a B. subtilis monolayer, and it would 
be very interesting to compare experimental rheological 
curves to our predictions for the rheology of active rolls 
and in-plane turbulent states. 
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Appendix A: Algorithms 

In our earlier work on these systems in ID P, the 
simulations for free and fixed boundary conditions were 
done using two entirely different algorithms. (This was 
for historical reasons.) After some initial benchmarking 
to check there were no systematic discrepancies between 
algorithms even in 2D, it has proved convenient to con- 
tinue with this division of labor, although in principle 
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FIG. 14: Top: Flow curve for C = 0.04, / = 0.002. Lower 
three panels: representative snapshots at a late time for 7 = 
0.00005 (second row), 7 = 0.001 (third row) and 7 = 0.002 
(bottom row). 
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FIG. 15: Effective viscosity (normalised by 77) versus time for 
an active gel with / = 0.002, ( = 0.0005, and fixed boundary 
conditions. The solid line refers to 7 = 0.0002, the dashed 
one to 7 = 0.001. The effective viscosities are smaller than 77 
in both cases, and nearly equal (demonstrating the existence 
of a linear regime). Note that in order to get an out-of-plane 
velocity field in steady state one needs to start e.g. with an 
out-of-plane director profile. 



either raethod should work equally well in both cases. 

Accordingly the work reported in this paper on free 
boundary conditions was done with a finite difference 
(FD) code developed from that previously used in [soj. 
That on fixed boundary conditions used instead a hybrid 
lattice Boltzmann (LB) code, in which FD for the order 
parameter is coupled to an LB description of momentum 
transport. This code was previously described in [31]. 

There are additional technical differences between 
these two sets of simulations. First, many of our LB 
simulations study square cells, Lx = 1.0, while our FD 
simulations study elongated cells, Lx = 4.0. Secondly, 
the FD work assumes a strictly 2D velocity field, so that 
the out-of-plane component Vz vanishes everywhere. In 
the LB runs, Vz{x^ y) can be nonzero although it vanishes 
at the walls via the no-slip condition there. This distinc- 
tion is only relevant when the system spontaneously ac- 
quires a flow along z. This may prove relevant in future 
work as discussed in Section IVTl 

We are interested in flows at low Reynolds number 
(Re) in which inertial effects are negligible. Within our 
FD simulations we therefore set p = directly, achiev- 
ing zero Re from the outset. In the LB approach, finite 
Re is inevitable since the algorithm exploits the physics 
of fluid inertia to evolve the dynamics: put differently, 
strictly zero Re requires strictly zero timestep. There- 
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fore, in the units defined above, the fluid density in the 
LB simulations is not zero but p = 2.3 X 10^. Despite 
this. Re values at the microscopic length scale / remain 
of order 0.1 or less in our simulations, small enough to 



maintain the inertialess character of the resulting flow 
(32I . (Even at the box scale. Re generally remains of or- 
der 1; moreover comparison with the FD results gives no 
suggestion that inertia becomes important here either.) 
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